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We study the properties of the subgap states in p-wave superfiuids, which occur at energies below 
the bulk gap and are localized inside the cores of vortices. We argue that their presence affects the 
topological protection of the zero modes. Transitions between the subgap states, including the zero 
modes and at energies much smaller than the gap, can alter the quantum states of the zero-modes. 
Consequently, qubits defined uniquely in terms of the zero-modes do not remain coherent, while 
compound qubits involving the zero-modes and the parity of the occupation number of the subgap 
states on each vortex are still well defined. In neutral superfiuids, it may be difficult to measure the 
parity of the subgap states. We propose to avoid this difficulty by working in the regime of small 
chemical potential fi, near the transition to a strongly paired phase, where the number of subgap 
states is reduced. We develop the theory to describe this regime of strong pairing interactions and 
we show how the subgap states are ultimately absorbed into the bulk gap. Since the bulk gap also 
vanishes as /i — > there is an optimum value fi c which maximises the combined gap. We propose 
cold atomic gases as candidate systems where the regime of strong interactions can be explored, and 
explicitly evaluate fi c in a Feshbach resonant 40 K gas. In particular, the parameter C2 parametrizing 
the strength of the resonance in such gases, sets the characteristic size of vortices, and the energy 
scale of the subgap states. 

PACS numbers: 74.20.Rp, 03.67.Lx, 03.75.Ss, 71.10.Pm 



I. INTRODUCTION 

Since the vision of a quantum computer based on the 
enigmatic degrees of freedom of topological phases of 
matter was set out P the search for physical realizations 
of topological phases has evolved as a leading topic of 
condensed matter physicsPThe physics of chiral p x + ip y 
BCS pairing^ represents a simple prototype system for 
topological order with prospective applications for inher- 
ently fault-tolerant topological quantum computing^ 

The p x + ip y paired phase is believed to occur in a 
number o f se ttings, including in the Al-phase of super- 
fluid 3 HeJ^ in the bulk of the layered perovskite ox- 
ide Sr2Ru04pEI as well as in two-dimensional samples 
of cold atomicP^ and polar molecular^ gases, and as 
a very closely related p-wave superconductor of compos- 
ite fermions in the v = 5/2 quantum Hall effectPSJEEU 
Physics similar to that of the topologically non-trivial 
Majorana modes in p x + ip y superconductors may also 
be induced by the proximity effect in interfaces between 
s-wave superconductors and topologic al insu lators,^ or 
related semiconductor heterostructuresJ^Hni 

In p-wave superfiuids, topologically non-trivial de- 
grees of freedom arise from the Majorana zero energy 
modes (ZEM) localized in vortices of the superfluid or- 
der parameter! 2 ^ The topological protection required for 
the ZEM to be used in quantum computing relies on the 
existence of an energy gap towards quasiparticle excita- 
tions. It is troubling therefore, that vortex cores in su- 
perfiuids feature eigenstat es o ccurring at energies much 
smaller than the bulk gapP^H 

In this paper, we discuss the implications of the pres- 
ence of such subgap states for implementations of topo- 
logical quantum computing (TQC) in a p x + ip y paired 



phase. We conclude that, while these states do not neces- 
sarily lead to decoherence of quantum information, they 
can complicate significantly the construction of any prac- 
tical scheme for TQC. In brief, at temperatures above the 
energy of the lowest subgap state ei , thermal excitations 
of the systems include processes which correspond essen- 
tially to a random flip of the qubit associated with the 
ZEM, represented by matrix-elements involving a single 
Majorana operator. However, as long as no excitations 
above the bulk gap are created, no decoherence may oc- 
cur, and information remains local to the vortex. There- 
fore, a compound qubit consisting of the ZEM as well as 
the complete set of subgap states is still well defined, with 
its state determined by the ZEM as well as the parity of 
the number of subgap excitations. 

The requirement to perform measurements of such 
compound qubits, however, may prevent implementa- 
tions of TQC in practice, particularly in neutral super- 
fluids where interferometry is not applicable. As a possi- 
ble solution to this dilemma, we suggest the use of spin- 
polarized atomic Fermi gases which may be driven to a 
strongly interacting regime of p x + ip y pairing in a con- 
trolled fashion (remaining in the weak-pairing phase^), 
by exploiting the physi cs of the BEC-BCS crossover^ in 
cold atomic gasesP2El This regime, where the bulk gap 
can be of the order of the Fermi energy was not consid- 
ered in detail in previ ous s tudies of vortex states in the 
Px + ip y paired phasesF^flMI I n particular, the existing 
scheme leading to an approximate analytical solution for 
the subgap stateiP does not apply. 

Following the need to elucidate the strongly paired 
regime of neutral superfiuids, in this paper we establish 
the theoretical framework for characterizing the physics 
of a quasi two-dimensional cold atomic gas in the BCS 
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regime above a p- wave Feshbach resonance. We show how 
to describe the system in a two-channel model, and de- 
duce the connection to the underlying parameters of the 
atomic gas, studying the example of 40 K, in particular. 
As a core result of this analysis, the size of the vortex 
cores is set by the parameter C2 that characterizes the 
strength of the resonance. We study the ensuing Bogoli- 
ubov de Gennes (BdG) equations in the limit of strong 
interactions and show how the subgap states merge with 
the bulk gap one by one until none is left. We then fo- 
cus particularly on the nature of the first excited subgap 
state, as opposed to a nu mber of previous studies that 
focused on the zero modes JMIEH22 Based on a numerical 
study of the BdG equations, we identify the regime with 
the best prospects for realizing TQC in p-wave superfluid 
cold atomic gases, i.e., the regime maximising the energy 
of the first subgap state t\ as a function of the detuning. 

The structure of this paper is as follows: in section \H\ 
we review the Bogoliubov de Gennes equations for the 
description of p-wave superfluids, identify the chemical 
potential as the single parameter of these equations in 
suitable rescaled units, and restate some known results 
about the subgap states in these units. Section [Til] pro- 
vides an in-depth discussion of how the subgap states 
influence the topological protection of the manifold of 
ZEM in a system with many vortices. We then proceed 
in section HVl to discuss the theoretical framework for the 
description of atomic p-wave superfluids in the regime of 
strong interactions, using a two-channel model. Solutions 
to the ensuing Bogoliubov de Gennes equations are given 
in Section |VJ mostly based on a numerical study with a 
strong focus on the properties of the first excited subgap 
state. Finally, in section |VD| we deduce concrete num- 
bers for experimental realizations of a p-wave superfluids 
based on a Feshbach resonance in potassium gases, before 



presenting further conclusions in section VI Details of 
how to extract the physical parameters for the Feshbach 
resonance in 40 K are given in Appendix |Aj Appendix [b] 
is devoted to a review of the approximate analytic cal- 
culation of the subgap states. Appendix [C] provides a 
discussion of the matrix elements between subgap states, 
and in Appendix [D] we review how to express the BdG 
equations on the sphere. 



II. MODEL 

Before entering the main discussion of this paper, let 
us introduce the notations of the formalism that we 
make use of below. The spectrum of spinless (fully spin- 
polarized) fcrmions whose p x + ip y pairing order param- 
eter is described by a gap function A(r) and gives rise to 
the Bogoliubov equations 



h — fl 7T 
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where, in a coordinate representation, the single particle 
kinetic term can be expressed as h = and the pair- 



ing term as 7r = |^/A(r) ( ^ - iJ^J ^/A(r). Note that 

in this definition of p-wave pairing, the gap function A 
has units of Energy x Length in contrast to s-wave super- 
conductors. Below, we adopt a dimensionless version of 
the Bogoliubov equations, obtained from ([I]) by rescaling 
the equations in terms of the characteristic length set by 
the gap function 



L = 



1 



mAo ' 



(2a) 



yielding the dimensionless length and energy-scales 



and E = EmL 



E 



(2b) 



Note the scale L is determined by the asymptotic value of 
the gap function in the bulk, A , measured far away from 
any vortices. We should note that while the gap func- 
tion Ao yields a characteristic energy scale Eq = mA§ of 
the problem, this scale is distinct from the bulk gap, as 
defined below, that represents a second useful reference 
energy. For the remainder of this paper, we use this di- 
mensionless formulation of the BdG equation 0, which 
is formally equivalent to setting the mass m — 1 and the 
value of the gap function in the bulk Aq = 1. This leaves 
p = /i/[mA§] as the single dimensionless parameter of the 
problem. In the remainder of this paper, symbols with a 
bar refer to the dimensionless versions, while bare sym- 
bols represent fully dimensional quantities. Wherever di- 
mensionless parameters appear in an equation, the other 
parameters are also dimensionless even if not explicitly 
indicated as such. 

Let us indicate a few known results in the dimension- 
less units. For example, the dispersion of the bulk quasi- 
particle excitations in the absence of any vortices now 
reads 



(3) 



The bulk gap A^ is set by the minimum of E^ which 
occurs at momentum at k = or at |fc| = \/2(p — 1) for 
p < 1 or /l > 1 respectively, with 



p. 



for fi < 1 



Ab VW^, for M > 1 



(4) 



Our study focuses on the spectrum of the p x + ip y su- 
perconductor in the presence of a radially symmetric 
vortex described by a winding of the order parame- 
ter according to A(r) = /i 2 (r)e 4K< ^, where h(r) —> 1 
at large r. Setting w(r) = exp[i(m + ^i)</>]M(r) and 
v(r) = exp[i(m — is =^)0]w(r), the Bogoliubov equations 
for general k read: 



-A + /I 



h 2 



d_ 

Or 




or 



v + hh'v = Eu 

u — hh'u = Ev 
(5) 
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where we introduced an abbreviation for the second order 
differential operator 

1 ~ 2*^ + ~2r~dr 2r 2 ' ^ 

In this paper, we focus on the cases of a vortex (n = 
1) and antivortex (k = —1). Analytical solutions 
to the equations ^ are known for several regimes of 
parameters P22 Most importantly, a zero mode exists for 
all vortices of odd vorticity, and is topologically pro- 
tect ed again st perturbations that do not destroy the bulk 
g a pf23M2ll j n a cylindrically symmetric vortex of vanish- 
ing size^and with vorticity k — 2n — 1, the wavefunction 
of the zero-mode takes the form 

with the (modified) Bessel functions J n (7 n ) [and the 
angular dependency involves and appropriate phase fac- 
tor as introduced above]. If a finite sized vortex core 
is considered, both the Bessel function and the exponen- 
tial localization are modified. The latter instead becomes 
exphj^V)*"]. 

It is als o kn own that there may exist additional sub- 
gap statesP^l with finite E < Ab known as the Caroli- 
deGennes-Matricon (CdGM) statesP^l While the energy 
of the zero-mode is protected, the energy of the CdGM 
states depends on the shape of the vortex. Their energy 
was calculated in the limit of p, 3> 1 for the closely re- 
lated case of a k = 1 vortex in the A phase of 3 HeJ^22l 
and found to be 

E m = -mui Q (8) 

for the subgap state with angular momentum m, as long 
as E <C Ab and with 

_ /-^exp^r^y)] 

"° ^dr^ V [-2f Q dr'h?{r')] ■ W 

Note the spacing of the subgap states is independent of 
the dimcnsionless ft, whereas in dimensional units loq ~ 
A B //j, <C A B P21 There are thus /x/A^ ~ y/JL such modes 
for large ft. Below, we consider the behaviour in the limit 
of small /2 and find that as p, is reduced, the subgap modes 
merge with the bulk one by one. For completeness, the 
perturbative solution for these eigenstates is included as 
Appendix [Bj Finally, we introduce t\ as a notation for 
the energy of the first subgap state, with ej — > luq in the 
limits of the approximation of large /2. 

III. ROLE OF THE SUBGAP STATES 

Given the presence of states below the gap energy, the 
question arises how these states affect the topological pro- 
tection of the Majorana zero-modes. The conventional 



view of the topological protection of a groundstate man- 
ifold requires the existence of an energy gap towards ex- 
cited states.^ At sufficiently low temperatures, the proba- 
bility of creating an excitation is then exponentially sup- 
pressed. If the vortices holding zero-modes are well sep- 
arated from each other, recombination would most likely 
occur in a way that restores the original groundstate con- 
figuration: due to the topological nature of the system, a 
quasiparticle excitation can only permanently change the 
state of the system if it braids around a second vortexP 
Excitations of the subgap states always remain localized 
to a single vortex. However, even if they cannot propa- 
gate information to a second vortex, they can make read- 
ing off the state of the qubit very complicated, as we see 
later. 

For p-wave superconductors in the BCS limit, the first 
subgap energy e± is typically much smaller than the bulk 
gap Ab (by a factor of Ag//i), and it may thus be im- 
practicable or even impossible to cool the system to tem- 
peratures below t\. For temperatures above ei, excita- 
tions of the subgap state(s) are likely. 

Transitions between the zero-modes and the subgap 
states require non-zero matrix elements connecting these 
states. It can be shown that matrix-elements for a scalar 
potential, created for instance by a passing phonon, are 
indeed non-zero by expanding a disorder potential V in 
the basis of Bogoliubov eigenstates. Therefore, transi- 
tions into the subgap states will occur at a finite rate in 
thermal equilibrium. In this case, the necessary energy 
for quantum jumps is then supplied by the heat bath, for 
example in the form of phonons in the ruthenatcs. An 
explicit calculation of the matrix elements involving the 
zero modes is shown in Appendix [C| For the purpose of 
this discussion, we only need to acknowledge the presence 
of processes of the form cj,7„ or j u c v , involving Majorana 

fermions 7„ and fermionic subgap states at vortex v. 

For any static external potential, the system has 
well-defined eigenstates with infinite lifetime. Time- 
dependent scalar potentials however, can provide the en- 
ergy to make transitions between eigenstates. The zero- 
mode has an interesting property in that its energy is 
protected against perturbations. However, its wavefunc- 
tion is deformed by a changing potential. Thus, the zero 
mode for a given external potential has a non-zero over- 
lap with excited states of an evolving potential at a later 
time. Non-adiabatic transitions can only be induced by 
perturbations, which occur sufficiently quickly on a time- 
scale set by the energy-scale for transitions, as the re- 
quired energy is supplied by the force resulting from a 
time-dependent potential. Therefore, the presence of the 
subgap states sets tighter limits on how stationary the 
scalar potential needs to be to conserve adiabaticity. In 
this context, we note that even a stationary disorder po- 
tential will act as a time-dependent perturbation driving 
non-adiabatic processes under braiding of the vortices. 

Having identified possible non-adiabatic processes 
among the zero modes and subgap states, the question 
remains whether these transitions cause decoherence in 
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the system. Let us first assume that we are able to 
momentarily cool the system to very low temperatures, 
even though non-adiabaticity cannot be avoided during 
braiding operations. Thus, we may start and end in 
the groundstatc. The consequences of temporary tran- 
sitions to the subgap states can be analyzed using the 
operator describing the braid which exchanges the po- 
sitions of two vortices 71 and 72, and which is given 
by P12 = -75(1 + 7172) (followed by a renaming of the 
vortices) 

For instance, excitation of the subgap state on vortex 
1, followed by an interchange of 1 and 2 and de-excitation 
of the vortex is equivalent to a simple exchange up a to 
sign: 

£272(1 + 7172)714 = -(1 + 7172) (10) 

Generally, the non-adiabatic processes affect braiding 
processes only up to a sign (or rather, up to a ran- 
dom abelian phase caused by the lack of knowledge of 
the energy as a function of time). More generally, any 
even number of transitions to subgap states will lead to 
an even number of Majorana operators being inserted - 
thus leaving the final state invariant. This is always the 
case if both the initial and the final state are within the 
groundstate manifold. However, should an odd number 
of subgap states be excited on a vortex, an error occurs 
and the final quantum state is altered [changing the en- 
tanglement properties of the wavefunction by modifying 
the relative phases of the components with empty and 
occupied core-states^j . 

Provided that all subgap states remain local to its 
vortex, and provided one can measure the number of 
fermions in the subgap states of all vortices, one can 
deduce how many Majorana operators must have been 
inserted in the initial state, and correct the final state 
accordingly. While such a procedure seems fundamen- 
tally possible, it would likely be extremely inconvenient 
in practice. 

It may be possible to find a measurement scheme which 
evaluates the state of the Majorana mode as well as the 
parity of the subgap states, which could be taken alto- 
gether as a compound qubit. Such qubits were considered 
independently in a recent paperPEI We now discuss the 
feasibility of measuring the atom number parity in the 
context of possible read-out schemes for quantum bits in 
p-wave superconductors / superfluids. 

As an example, let us consider a proposal to detect 
the state of qubits in atomic p-wave superfluids via 
spectroscopy.^ This scheme relies on the possibility to 
detect the presence or absence of a single unpaired atom 
after fusing two vortices forming a qubit. Detecting a 
single unpaired atom is possible as the coupling of the 
atom's internal states to a Raman pulse depends sensi- 
tively on the detuning of the incident radiation. How- 
ever, the required detuning differs by twice the bulk gap 
between paired and unpaired atoms, enabling one to ad- 
dress only the latter. Given the presence of the subgap 
states, it is now possible that there are multiple unpaired 



atoms in a single vortex. Consider, for example, a situ- 
ation where an even number of subgap states have been 
excited in the two vortices of the qubit. Then, according 
to the proposal under consideration,^ we would detect 
the presence of unpaired atoms; however, the state of the 
qubit corresponds to an empty |0) state. To deduce the 
state of the qubit, it is therefore required to detect the 
number of unpaired fermions (which equals the sum of 
occupation numbers of all subgap states) to be able to 
deduce its the qubit state from the parity. This requires 
more accurate detection techniques. 

Interferometry appears as a suitable measurement 
method, as it is sensitive only to the parity of the num- 
ber of atoms. In charged superfluids, interferometry is 
expected to workP provided that transport along the 
edge can be described in terms of the motion of quasi- 
particles with no internal degrees of freedom. However, 
in compressible neutral p-wave superfluids, interferom- 
etry is certainly not applicable at all as the compress- 
ibility of the system entails particle number fluctuations 
that wash out the interferometric signal when braiding 
quasiparticlesP^ For incompressible states such as the 
Moore-Read phase found in the v = 5/2 quantum Hall 
effect, interferometry is considered as a leading scheme 
for read-out P we should also note in passing that there 
are no subgap states in the spectrum for FQHE states, as 
the size of the vortex is of the order of the interparticle 
spacing. 

As the presence of the subgap states gives additional 
structure to the quantum bits, we anticipate that the 
read-out protocols will be more complex, and thus be 
prone to errors in the read-out process, particularly for 
neutral superfluids. Therefore, the presence of low-lying 
subgap states very likely prevents the practical use of the 
quantum register in these systems. 



IV. P-WAVE SUPERFLUIDS IN THE LIMIT OF 
STRONG INTERACTIONS 

In the previous section, we argued that the presence of 
low lying subgap states complicates the use of p-wave su- 
perfluids for TQC. From a technological point of view it 
seems imperative to evade the problem of having to keep 
track of the occupation numbers of all subgap states. 
In p x + ip y superfluids of cold atomic gases, this can 
be achieved by driving the system into the regime of 
strong interactions where the number of subgap states 
N s is expected to be small. For Jx 1, the semiclassical 
approximation^ yields N s cx yfp,. The core task of this 
paper consists in establishing the nature of the subgap 
states for p, small. 

To do so, we need to solve the BdG equations for a 
vortex when p, is small. This can only be done numeri- 
cally. The appropriate equations depend on the shape of 
the condensate h(r), which vanishes linearly in r near the 
core for singly quantized vortices. At the same time, as 
r increases h(r) approaches 1 on a characteristic length- 
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scale called the coherence length of the condensate, £. 
We need to know £ to fully define the BdG equations. 

A known way to create p-wave superconductors with 
p, which is not very large is by employing cold atomic 
gases with p-wave Feshbach resonances. Current exper- 
imental realizations of atomic gases with p-wave Fes- 
hbach resona nces a re not stableJ^HSl due to 3-body 
recombination PHHD However, there are p roposals to en- 
hance stability using optical lattices P^IMI j n |- ne nex ^ su )- ) _ 
section we study such systems in order to determine its 
parameters, including the scale of the gap function, the 
chemical potential and the coherence length of the con- 
densate. Additionally, wc also consider p-wave supcrflu- 



ids of dipolar molecules in section IV B 



the formation of the molecules, typically about 1 — 3nm. 
Therefore the formation of molecules remains a 3D pro- 
cess. 

The detuning S, the parameter which can be controlled 
experimentally by varying the magnetic field, allows to 
vary the strength of pairing between the fermions. 6 is 
simply related to the "bare detuning" eoj the parame- 
ter which appears explicitly in the Hamiltonian, by the 
relation 



6 = 



eo — const 

1 + C 2 ' 



(12) 



where const is an irrelevant constant, and ci is a param- 
eter which characterizes the strength of the resonance, 



A. Atomic Fermi gases 

In the context of cold atomic gases, p-wave pairing 
is studied via the two channel model that considers 
free fermions and their bosonic bound states as separate 
entities. 12 Its Hamiltonian reads 



p q,« 

p,q,a 



+p 



Am J Qq Qq 



at 



h. 



(11) 



This model describes a gas of fermions with creation 
and annihilation operators a\ a which can form bosonic 
molecules with the angular momentum 1 (hence the 
bosonic molecules are described by the creation and anni- 
hilation operators w a , b a with the vector index v). Once 
the bosonic molecules Bose-condense, the fermions form 
a p-wave superconductor. For homogeneous p-wave su- 
perfluids, it was found t hat t he p x + ip y paired phase 
is always the groundstateP^ T n the experimentally rel- 
evant case of gas with a dipolar anisotropy, the triplet 
of m = +1,0,-1 states is split, and the m — state 
becomes the groundstate, and as a result p^-pairing is 
present in the phase diagram. For weak anisotropies, the 
chiral state remains the low-temperate phase over a large 
range of detunings. Even in the case of strong splitting, 
the chiral p x + ip y paired phase remains present: Fesh- 
bach resonances for the m = and m = ±1 channels are 
then also well separated, and the chiral phase is observed 
near the latter oneP^l 

We observe that in real experiments, the cold atomic 
gas would be confined to two dimensions up to a "pan- 
cake" of width I of the order of the wavelength of visible 
light, or 500nm. The momenta in Eq. (Ill are chosen 



appropriately to reflect such a geometry. We call this 
setup a quasi-2D gas (unlike a purely 2D setup where all 
motion is completely confined to 2D geometry). 

The confinement length £, while much smaller than 
interparticle separation, thus leading to a truly 2D su- 
perconductor, is still much larger than the molecular size 
R e , which is set by the range of the forces responsible for 



to 2 to g 
C2 = ^ 9K =Z^R P 



(13) 



The parameter A = l/R e is hidden in the two-channel 
model, appearing as an upper cutoff for all the sums over 
momenta. 

The solution to the two-channel model can be de- 
scribed in the following way. First of all, the bosons are 
always Bose condensed, or 



(14) 



VB a 



This leads to fermions forming a p-wave superconductor. 
Second, the vector structure of B a defines the order pa- 
rameter, and for a,p x +ip y pairing, energetically favorable 
within this model, 



B x — —iBy — B. 



(15) 



Given the density of particles n and the corresponding 
Fermi energy (understood as the Fermi energy of a free 
Fermi gas at this density), we can identify three regimes. 
If S > 2ep, then the density of bosons is exponentially 

smalP^lin the parameter exp |~ '*~^,' :F «s}i where S is de- 



fined below in Eq. (30 1. Since the density of bosons 



n b = 2B 2 



(16) 



is responsible for the superconducting pairing in this 
problem via 



A = 2gB 



(17) 



the regime of exponentially small rib corresponds to the 
conventional BCS superconductors in which the transi- 
tion temperature is a small fraction of the Fermi energy. 
The chemical potential in this regime is fj, = ep ■ 

Next, if < 8 < 2ep, a finite fraction of fermions 
converts into bosons, rib is now of the order of the initial 
density of fermions in this problem, and the superconduc- 
tor which forms under these conditions has a transition 
temperature which is a substantial fraction of the Fermi 
energy (and growing as 6 is decreased). At the same time, 
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the chemical potential of the fermions p is approximately 
equal to (5/2, 



(18) 



where "approximately" means up to terms of the order 
of g 2 ). Thus in the terminology of Ref. [3] this is still a 
weakly paired superconductor (that is, a superconductor 
with p > which has large Cooper pairs). We refer to 
this regime as strongly interacting. 

Finally, when S becomes negative, the chemical poten- 
tial changes sign. Now the density of fermions is exponen- 
tially small, while most particles are bosonic molecules 
which are Bose condensed. The transition temperature 
of such a superconductor is a certain finite fraction of 
the Fermi energy. In the terminology of Ref. [3| this is 
now a strongly paired superconductor, separated from a 
weakly paired superconductor by a quantum phase tran- 
sition which occurs at \i = or 5 close to 0. 

Returning to the discussion of possible experimental 
realizations of two dimensional gases with cold atoms, 
the two-channel model describes a two dimensional su- 
perconductor which, at fi > (implying positive 5), is 
the topological state of matter of interest to us, here. 
A typical experiment would be conducted in the regime 
where < S < 2cf to maximize the transition tempera- 
ture of the superconductor, and as we show later, to re- 
move the subgap states for a particular S. In this regime, 
< p < ep, unlike the conventional superconductors 
where ju, = ep. The principal goal of this paper is to 
understand what happens to the vortex subgap states as 
the chemical potential becomes smaller than ep. 

We now use this scenario of the quasi-2D two-channel 
model to estimate the coherence length of the condensate 
£, which we need to be able to estimate the typical size 
of a vortex. 

To do that, we integrate out the fermions and concen- 
trate on the effective action of bosons. This was done in 
Ref. fT2l with the result, in the regime where \i < ep, 



S 



dVdt 



Cg c 2 m 



4m 



M)Ml + c 2 )+ 
(19) 



Here C is a dimensionless constant whose precise value 
is not important for our purposes here. We note that 
the calculations in Ref. are done in 3D, while we work 
in quasi-2D. Howev er, t he contributions from integrating 
out the fermions in ( 19 ) are all proportional to c 2 ~ g 2 A. 



In other words, they come from the momenta of the order 
of A ^ l/R e 3> 1/t where I, the width of the condensate, 
is much larger than R e , the range of the interactions. So 
Eq. |l9| is valid in quasi-2D, as well as in 3D. 

The coherence length £ can be extracted by comparing 
the kinetic and quartic terms of Eq. (19). We find 



l + c 2 
m£ 2 



g c 2 mB , 



(20) 



where we replaced J2 a bab a by B 2 in the spirit of Gross- 
Pitaevskii equation. This gives 



(-2 



gBmy/c2 



(21) 



Throughout this paper, however, we are interested in 
the dimensionless £, expressed in units provided by A 
[the units of length given by l/(Aom), see Eq. ^2|]. In 
turn, the dimensionless coherence length can be found as 



C-tA \A + C 2 o 

? = ?^o™ ~ —p, —mgB 

gBraJc2 



1 + C 2 

C2 



(22) 



Thus we find that if c 2 is large, the coherence length in 
our units is close to 1. If c 2 is small, the coherence length 
can be larger than 1. 

In the p-wave Feshbach resonance in 40 K, c 2 can be 
shown to be around 14.4 (see Appendix [A]), and the co- 
herence length is close to 1. However we do not know the 
value of C2 for other p-wave Feshbach resonances, so we 
cannot make any assumptions about it beyond Eq. (22 1. 



We can further use (11 1 to calculate how the bulk gap 



Ao and the chemical potential /i depend on the detuning 
6. To do this, we rely on the following arguments. 

If 5 > 2ef, then Ao is exponentially small. At the same 
time [x — ep. When expressed in units of mA 2 ,/^ 2 , as we 
do throughout this paper, p, 1. This is the regime of 
conventional superconductors. 

When S is lowered below 2ep, then A quickly starts 
to increase. Let us determine Ao as a function of 5 in 
this regime. As before, we work in the quasi-2D regime, 
where the p-wave gas is confined to a pancake of width 
£, such that I is much smaller than the average parti- 
cle separation. We write down the particle conservation 
condition 



E 



p 

2 in 



PL 

2m 



AoV 



2 • 2B 2 V = Ni 



total ■ 



(23) 



Here, B is the condensate density originally defined in 
(14), and Aq = 2gB. The summation over p reflects the 



quasi-2D geometry of the p-wave gas and may not be 
straightforward to convert into an integral, as is usually 
done in these cases. However, we observe that the sum- 
mation over p is actually divergent at large p~A>l/f. 
This allows us to capture this divergence by introducing 
a 3D integral J2 P ~^ V I d 3 p/(2ir) 3 . This was already 
done in Ref. [T2J with the result 
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In this equation, the summation over momenta p is now 
restricted by p <C 1/2, while the contribution of the lead- 
ing divergence at larger momenta can be absorbed into 
the B 2 term. The remaining summation over p is purely 
two-dimensional. We can evaluate it by approximating 
the expression in the square brackets as a Fermi-Dirac 
step. This yields 



V m/i 



4B 2 V(l + c 2 ) = Ni 



total ■ 



(25) 



Note the volume V is a three-dimensional quantity, which 
yields the corresponding 2D volume as V/£. To express 
our final results, let us write also 



5 
2ep 



J/_ 



(26) 



which varies from to 1 and measures detuning 8/2 in 
the units of Fermi energy. In terms of this parameter, 



B = 



'ej7m(l — x) 



8tt£(1 + c 2 ) 
Finally, this gives Aq = 2gB as 



A = g\ 



'ei?m(l — x] 
27r^(l + c 2 )' 



(27) 



(28) 



For the purposes of this paper, we would like to compute 
fi = 5/2 = xep in the units of mAg. This gives the 
dimensionless chemical potential 

fi= -A^i =27r(l + c 2 )^-^. (29) 
mAg m z g z 1 — x 

For future reference, we name the pref actor in this equa- 
tion which sets the overall scale of the gap and chemical 
potential 



S = 2tt{\ + c 2 ) 



m 2 g 2 



(30) 



We can see that as the detuning is decreased past 2ep, x is 
varied from 1 to 0, the dimensionless /2 indeed varies from 
very large values [( 29 ) predicts infinity at x = 1, although 



this is an approximation artifact; at large detuning we are 
in the regime of conventional superconductor with a very 
large but finite p] all the way down to 0. 



B. Polar molecular Fermi gases 

For diatomic molecular gases in two dimensions, 
attractive interactions may be generated by dress- 
ing the molecules with circularly polarized microwave 
radiation.^ The result is a dipole moment which is ro- 
tating in the plane of the 2D gas, and the interaction 
averaged over the angle of rotation yields a net attrac- 
tive long-range potential V{r) = -d 2 ff /(2r 3 ). (We follow 



the notations of Ref . IT3"1 ) The strength of this interaction 
depends on the field strength for the incident microwave 
interaction, as well as on the permanent dipole moment 
of the molecules. It is characterized by a lengthscale 
r* = Md 2 cS /2h 2 , that can be of the order of r* ~ 200nm 
for realistic experimental parameters in typical 7 Li 40 K 
molecules. The dimensionless strength of the interaction 
kpr* , can therefore be of order one in gases of a fairly 
low density.^ 

For superfluids with dipolar interactions, theory has 
been developed only at the level of a BCS mean-field 
description! 1 ^ I n this framework, the critical temperature 
and bulk order parameter Aq are obtained to be 



(31) 



For large kpr*, we expect that a significant fraction of 
the fermions will be paired, leading to a reduction of the 
chemical potential for fermions. However, to quantify 
this effect more theory would have to be developed to 
study the specific case of dipolar interactions. 

Similar to the reasoning leading to Eq. (20), the 



coherence-length £ is obtained by balancing the kinetic 
and interaction terms in the underlying pairing Hamilto- 
nian. The result can be expressed in terms of the Fermi- 
velocity and the critical temperature 



\pe 



(32) 



and when stated in our dimensionless units, this equates 
to a coherence length of order one 



£ = £ m Ao = Xpmep /k F — !• 



(33) 



Note that generically, £ ~ 1 for any superfluid in the BCS 
limit, as we have used no specific features of the dipolar 
interaction to derive this relation. 



V. SOLUTIONS OF THE BDG EQUATIONS 

The two-channel model discussed in the previous sec- 
tion provides a framework to discuss the physics of 
molecule formation in strongly interacting superfluids. 
Formally however, in strongly interacting p-wave super- 
fluids with large C2, one can integrate out the bosons 
in this problem and reduce the two-channel model to the 
one-channel mode l desc ribed by the usual Bogoliubov de- 
Gennes equations! 12 * 43 ! In this section, we discuss the so- 
lutions of this equation, focusing on the properties of the 
subgap states for p, < 1. We first point out some gen- 
eral features of the analytic solution, and then deploy the 
formulatio n of t he Bogoliubov de Gennea^ equations on 
the spherep^ll as we ll as its numerical solutions. 



A. Asymptotic solution of the BdG equations 

As an approach to discussing the solutions of the BdG 
equations for a cylindrical vortex ^ , let us first analyze 
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FIG. 1: (color online) Inverse localization length £ 1 as a 
function of the energy of a core state. 



the asymptotic behaviour in the limit of r — >■ oo. The 
solutions are known to be strongly oscillatory functions, 
so the derivatives of the Bogoliubov functions are of the 
order of the functions themselves. However, terms in 1/r 
can be neglected, as well as h'(r) = far outside the 
vortex. Thus, ^ reduces to 
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(34) 



Using the Ansatz of (u v) T — (A B) T exp{i7r}, and 
solving the characteristic equation for 7 at a given value 
of the energy yields four possible solutions of the form 



7±,± 



± 



2(p - 1)± 2i 



E 2 



(35) 



First, note that in the limit of E — > 0, p 1, these sim- 
plify to 7±,± ~ ±v / 2/2T*- Among the two solutions which 
are finite as r — > 00, we recover the behaviour of the 
modes Q, namely oscillations with wavenumber ~ kp 
and an exponential decay with a characteristic length- 
scale of localization Iq = 1. In the general case, let us 
decompose 7 = k(E) + Uq 1 (E), which, for p > 1 can be 
written as 



k(E) = ±\ 2dp 2 - E 2 cos 



1 

- arctan 
2 p-1 



E 2 



E 2 t 



- arctan 
2 p 



E 2 



(36) 



According to these equations, the wavenumber of the vor- 
tex state depends only very weakly on the energy, varying 
between fc(0) = y/2p- 1 and k(A B ) = ^J2p - 2. This ef- 
fect is significant only if p is small. On the other hand, 
the localization length is one at zero energy ^o(O) = 1 and 
diverges as E — » . The dependency of the localization 
length on the energy of the subgap state is illustrated in 
Fig. [T] for several p. For p large, the localization length 
goes as t = [1 - (E /A B ) 2 }^ 2 . 



B. BdG equations on the sphere 

To study the physics of vortices in a finite size system, 
a convenient choice is to place a vortex-antivortex pair on 
the surface of a sphere. The BdG equations for this con- 
figuration were recently described by Kraus et aiP^ They 
are obtained by expanding the p-wave pairing function on 
the spher e A (1 7, Of) in terms of the spherical monopole 
harmonic d 47 * 48 * Y^ 9 m , which also serve for the expansion of 
the Bogoliubov eigenfunctions u, v. For convenience, we 
provide the resulting equations in dimensionless units in 
Appendix [DJ The relevant dimensionless parameters are 
given by the radius of the sphere R, and the size of the 
vortex core £; both are measured in units of the length 
scale L from Eq. [2j The dimensionless chemical potential 
p is set by the Fermi angular momentum lp (number of 
shells filled) , and also depends on the radius of the sphere 



If(If + 1) 
2R 2 ' 



(37) 



C. Numerical results 

In this section, we discuss the spectrum of Bogoliubov 
quasiparticle excitations on the sphere with a vortex- 
antivortex pair, which is obtained by solving the eigen- 
value problem in its matrix form (Dl ) with standard lin- 
ear algebra packages. 



1. Global Spectrum 

In Fig. [2j we display the energies of eigenstates found 
below the bulk gap As as a function of the chemical 
potential, while areas above the bulk gap are shaded blue. 
These spectra were calculated for a sphere of radius R = 
40, and include values of the Fermi angular momentum 
If = !)•••; ' w ith a cut-off for the equation set at 
( mal = 200. The bulk spectrum opens like a cone near 
p = 1 and crosses over into a square root behaviour, 
following Eq. Q. 

The energy of the subgap states found to depend 
weakly on the chemical potential. For the Caroli-de- 
Gennes-Matricon (CdGM) statesf^the prediction is that 
the eigenstates have strictly no dependency on p, as well 
as a linear dispersion in angular momentum [see Eq. ([8])]. 
Our results confirm that this is true for large enough 
chemical potential p and vortex size £. In the bottom 
panel of Fig. [2] the energy of the first subgap state ei 
remains roughly constant from p = 2 down to the point 
where it is absorbed into the bulk. For vortices with 
smaller cores, as shown in the top and centre panel of 
that figure, there is some 'bending' of the subgap states: 
the energy of the subgap states goes towards a constant 
only for large p, while at smaller chemical potential the 
absolute value of their energy decreases, smoothing into 
the cone of the propagating Bogoliubov quasiparticles. 
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FIG. 2: (color online) Spectra of a p-wave superconductor 
with a vortex-antivortex pair at the antipodes of a sphere 
with radius R = 40 as a function of the chemical potential 
p. From top to bottom, the figure shows the spectrum for a 
small vortex with £ = 1, a vortex with £ = 3 and a vortex with 
£ = 10. Blue shaded areas indicate the presence of excitations 
above the gap. The subgap states are marked by crosses and 
include excitations below the gap up to a angular momenta 
with \(L Z ) — ~| < 10. We use the dimensionless units of 
Eq. The vertical line at p « 0.6 situates the cut through 
the spectrum, which is shown in Fig. [3] alongside p w 2. 



However, the energy where the subgap state is absorbed 
into the bulk spectrum is still of the same order as the 
asymptotic value. 

In addition to this slight p dependency, another feature 
is apparent in the spectra: while low-lying subgap states 
are evenly spaced, the density of states exhibits a step 



increase at a value below the gap, signalling the presence 
of additional subgap states. 

Another phenomenon, best visible in the top panel 
with £ = 1, is the splitting of the subgap states. Each 
of the subgap states occurs as a doublet consisting of 
one state localized in cither pair of the vortex and an- 
tivortex present on the sphere. While the splitting of the 
zero-modes results from the hybridization of two modes 
found precisely at zero energy, the splitting of the modes 
at non-zero energy is predominantly of a different na- 
ture: as we show in more detail below, the wavefunc- 
tions in the vortex-core are distinct for the vortex and 
antivortex-state for E ^ (the splitting of the doublet 
saddling E = is not visible in this figure). 



2. Spectra at fixed fi 

The dispersion of the spectra as a function of the an- 
gular momentum to — (L z ) relative to the symmetry 
axis of the vortex cores reveals some additional insights. 
For simplicity, the discussion focuses on the states with 
E > 0, keeping in mind that states at negative energy 
are related by virtue of the symmetry of the Hamilto- 
nian E(m) = —E(l — m). We display the dispersion 
at two distinct values of the chemical potential. Figure 
[3] shows the energy as a function of to for p « 0.6 (left 
column) and p ~ 2 (right column) respectively, in the ge- 
ometries already used for Fig. [2] (corresponding to values 
of l F = f and l F = ±p). 

Let us first discuss the case of p — 0.6. Irrespective 
of the size of the vortex core, we find that the bulk gap 
Ab (solid red lines) conforms well with the case of a ho- 
mogeneous order parameter without vortices, as given by 
Eq. Q. By contrast the estimate for the energy of the 
first subgap state (black dotted li nes) , obtained from ^ 
using the same radial profile as in (D2 1, is much less accu- 



rate, in particular for the small vortex cores. For £ = 1, 
the predicted subgap energy wo(£ — 1) ~ 0.814 is roughly 
1.5- fold larger than the actual eigenstates l\. The chem- 
ical potential here was chosen such that the first subgap 
state almost merges with the bulk spectrum. The cor- 
responding pair has a very large splitting of S± = 0.042, 
almost 8% of their median eigenvalue e\. The splitting 
of the states is analysed in more detail, below. With in- 
creasing size of the vortex core, the previous estimate of 
ei Riwo from Eq. ([9| is increasingly accurate. However, 
in this regime of small p, the subgap states have a dis- 
persion which is sub-linear, i.e., the mode of the CdGM 
states bends to asymptote the bulk gap. 

In the bottom panel of Figure [3j it is apparent that 
there can be multiple subgap states at a given angular 
momentum m. Like the CdGM states, these modes con- 
sist of pairs of eigenstates at each value of to, as ex- 
pected for localized states in the presence of two vortex 
cores. We refer to these states by an additional quan- 
tum number n, where n = refers to the CdGM branch. 
In addition, we denote the two degenerate states of each 
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FIG. 3: (color online) Spectra as a function of the angular momentum m = (L z ) for a p-wave superconductor with a vortex- 
antivortex pair at the antipodes of a sphere with radius R — 40 at fixed chemical potential fi. The left column shows fi w 0.6 
(according to the yellow cuts in Fig. 2j, and the right column shows fi — 2. Vertically, we show configurations for a small (top, 
£ = 1), medium (centre, £ = 3) and large (bottom, £ = 10) vortex respectively, as in Fig. [5] The expected bulk gap As is 
shown as a solid red line, the energy ujo according to |9| as dotted black lines. The subgap states can be classified as two types: 
the anomalous branch of the CdGM states are approximately linearly (for small m) dispersing near E — 0, with a negative 
slope. Additional subgap modes with radial quantum number n / occur in large vortices and have with a finite minimum 
(maximum) of \E\ ~ O(Ab) at/near m — \. The number n ^ modes increases with p,. 



twofold degenerate branch by +/— . The dispersion of the 
)i ^ modes has a minimum in \E\ at or near m = i. 
On first sight there appears to be a symmetry relating 
eigenstates of E{m) o E(l — to), but the spectrum is 
slightly skewed such that E(i + 8m) > 2?(~ — 5m), given 
(6m = 1, 2, . . . > 0). This is unlike the case of s-wave su- 
perconductors where the n =^ modes are characterized 
by a symmetry of the spectrum in E(m) = E(—m) P^For 
large |m|, the modes asymptote towards the bulk gap. 



Near their minimum, the dispersion of the n ^ modes 
is roughly quadratic, leading to the markedly higher den- 
sity of states as compared to the n = CdGM states, and 
thus explaining the jump in the density of states that we 
pointed out in Fig. [2] 

We are not aware of other numerical studies that have 
analyzed the n ^ subgap states, as in typical type II 
superconductors, the vortex core is too small for such ad- 
ditional bound states to exist Within BCS theory, we 
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found above in section [IV B| that the coherence length of 
the superfluid is generically near one. In atomic p-wave 
superfiuids, the existence of the dimensionless parameter 
C2 allows one to scale the coherence length, which may 
become large for small C2 according to Eq. (22). There- 
fore, discuss their properties in more detail in section 
[VCll below. 

Due to the asymmetry of the n^O subgap modes be- 
tween positive and negative angular momenta, they may 
influence the dynamics of the vortices in superfiuids. In 
particular, the Magnus-force acting on vortices that are 
moving through the system may acquire corrections re- 
sulting from the influence of the subgap states EESISHES] 
although it was debated whether such influences should 
vanish as the Magnus force can be derived from a Berry 
phased In s-wave superconductors, the CdGM states 
would be singled out as the only mode contributing to 
corrections to the Magnus force, as these corrections van- 
ish if ^E n (k) is oddP^In p-wave superfiuids, J^E n (k) 
does not have this symmetry and thus, n^O modes may 
contribute. It is difficult to estimate the magnitude of the 
additional corrections: while the asymmetry of J^E n (k) 
is small, the density of states contributing to corrections 
is much higher than for the CdGM mode. 




FIG. 4: (color online) The wavefunctions u(r), v(r) for the 
first excited state of the CdGM branch with n — and 
m = 3/2, for a sphere of radius R — 10 and chemical po- 
tential p = 8. The eigenvalues associated with the upper and 
lower eigenvalue are split and amount to eo+ = —0.441639 
and eo- = —0.443256, respectively. The upper panel shows u 
and v superposed with a fit that shows excellent agreement 
for the behaviour in the vortex core. The lower panel dis- 
plays the wavefunctions on a logarithmic scale and indicates 
the expected exponential decay for comparison. 



3. Splitting of the CdGM states 

In section |V C 1| above, we noted that the splitting of 
the CdGM states in the spectra for the vortex-antivortex 
pair is much larger at E ^ than for the Majorana 
zero-modes. For the latter, this splitting is interpreted 
as arising from a tunnelling term that hybridizes the two 
degenerate modes. 

Contrary to the case of the zero modes, subgap states 
at m 7^ 5 ( or E 7^ 0) are not precisely degenerate even 
for a well separated vortex-antivortex pair. This may be 
suprising at first, as within the perturbative solutiorPof 
the BdG equations, the energy of the subgap states does 
not depend on the vorticity. 

However, the BdG equations for different vorticity 
K are distinct: as Eq. ^ shows, the equations dif- 
fer in the index of the differential operator 2?;, result- 
ing in different short distance behaviour via the terms 
[m ± (k — l)/2]/r 2 . Indeed, the perturbative solution 
yields wavefunctions proportional to the Bessel functions 
u m , K (r) ~ J m+ n=±{r) and u m ,«(r) ~ J m _^i(r) at small 
r. According to this solution, u and v have the same short 
distance behaviour behaviour in a vortex, while they are 
proportional to two distinct Bessel-functions in the an- 
tivortex, with an index differing by two. (Only for the 
zero-mode this offset is trivial as it pairs u ~ J\ and 
v ~ J—x, which are equal up to a sign.) It then becomes 
obvious that the energy of the eigenstates depends on the 
detailed interaction between the shape of wavefunctions 
and the order parameter within the vortex-core that is 
neglected in the perturbative solution. For illustration, 
Fig. [4] displays the eigenfunctions for the first subgap 




2 4 6 8 10 12 1 10 

R 5 

FIG. 5: (color online) Left panel: The splitting of the zero 
modes and first CdGM subgap states at m = 3/2 is shown as 
a function of the radius R of the sphere geometry for several 
values of p. The splitting of the zero modes follows the am- 
plitude of the envelope of the wavefunction at R/2 (magenta 
dashed). For small p, some oscillations are seen, as displayed 
separately below in Fig. [6] At large enough vortex separation 
r = ttR, the splitting of the m = 3/2 states is given by the 
shape of the vortex core and independent of r. At small r 
it is increased by the hybridization of the modes. Top right 
panel: splitting of the subgap states 5E at m = 3/2 as a func- 
tion of p, found to be inversely proportional. Bottom right 
panel: the same splitting SE(m = 3/2) as a function of the 
vortex size. The dashed lines show a fit with a third order 
polynomial in 
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state, showing clearly the distinction between the cases 
of a vortex and antivortex. 

A detailed analysis of this situation is presented in 
Fig. [5j where we focus in particular on the splitting of 
the zero mode, as well as the first subgap state. In the 
left panel, the magnitude of these splittings is plotted 
as a function of the radius of the underlying sphere R, 
which sets the distance of the vortex-antivortex pair. The 
splitting of the zero mode is described very well by the 
square of the envelope of the wavefunctions at half the 
sphere radius exp[— 2 J^ 2 h 2 (r)dr] w exp[— 2.R], for low- 
lying subgap states. We discuss some deviations observed 
at small chemical potential in more detail, below. For 
the first excited state, the splitting falls onto this curve 
only if the radius is very small. At large separation, it 
saturates as a constant value and is entirely set by the 
physics inside the vortex core, as argued in the previous 
paragraph. The two panels on the right-hand side indi- 
cate the magnitude of the splitting for the first subgap 
state as a function of either the chemical potential p and 
the vortex size £. The functional dependency goes as 
5E(m = 3/2) ~ p- 1 - Its relationship to the size of the 
vortex core is less clear, but it can be fit well by a third 
order polynomial in the inverse core size, in accordance 
with the notion that the shape of the vortex core strongly 
influences this energy scale. 

The amplitude for tunnelling processes between vor- 
tices is set by the overlap of the respective wavefunctions. 
As the eigenstates are strongly oscillatory, the tunnelling 
amplitude as a function of the vortex separation is not 
merely set by the exponential envelope, but it addition- 
ally changes sign periodically in the separation betwe en 
vortices, as recently discussed in the literature! 28 * 55 ! It 
may be surprising at first that no such oscillations are 
seen here. However, the configuration we study is very 
special in that the circumference of the sphere is close an 
integer multiple of the Fermi-wavelength for large lp ac- 
cording to (37). Deviations from this situation are thus 



found only at small values of the Fermi angular momen- 
tum. Indeed, we find that oscillations can be seen for 
small lp, as shown in Fig. [6] Given that the chemical 
potential can take only the discrete set of values (37 1, 



we plot the splitting at constant lp, rather than constant 
p. Note that the oscillations wc see occur at a lower 
frequency for larger chemical potential, even though the 
wavelength of oscillations for the wavefunctions of these 
states increases, in line with our explanation in terms 
of multiple Fermi-wavelength filling the circumference 
of the sphere. The oscillating behaviour of the subgap 
states was recently explored numerically in the plane 
geometryP2l 



4- Subgap states: branches with n/0 

In large vortex cores, multiple branches of subgap 
states occur, which can be interpreted as states with 
a different radial quantum number n (the CdGM being 
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FIG. 6: (color online) Splitting of the zero modes 5Eo at m = 
| for small values of the Fermi angular momentum ranging 
from lp = 12 (leftmost panel) to If = 24 (rightmost panel), 
plotted as a function of the sphere radius R. See main text 
for a discussion. 
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FIG. 7: (color online) (a) The energy £1/2,1 of the eigen- 
value at angular momentum m = s for branches of subgap 
states with radial quantum number n = 1 (solid lines) and 
n — 2 (dashed lines), for different chemical potentials p,. This 
data was collected for a sphere with radius R = 8. When 
expressed in units of the bulk gap As (as shown), the de- 
pendency E/Ab approximately collapse onto a single curve 
for large p,. (b) Showing the dependency of £1/2,1/ As on 
the chemical potential for the n = 1 branches of the subgap 
states: the additional branches are pushed into the continuum 
for small fl. 



n = 0). Here, we analyse some properties of the eigen- 
states at angular momentum m — \ at different values 
of £ and p, with eigenvalues denoted as E m n . 

As displayed in Fig. [7J our numerical results show that 
additional modes occur for vortices with size £ > 1. The 
energy of the n > 1 branches can be expressed in fractions 
of the bulk gap: to a good approximation, the eigenvalues 
do not depend on the chemical potential (except for small 
p), and Ei n (^,p) = f n (£)A B (fi). 

We now discuss features of the wavefunctions of the 
n > modes. Firstly, as these states occur at energies 
which are large fractions of the bulk gap, they are less 
strongly localized than the low-lying subgap states. An 
example set of wavefunctions for a vortex pair with £ = 5 
at p = 8 on a sphere with R — 12 is shown in Fig. |8j 
where the upper panel gives a comparison of the n = 1 
state (occurring at E^mi/Ab ~ 0.6) and the n = 
zero mode in logarithmic scale, showing the excited state 
with a localization length of about £q — 1.35, compared 
to £q = 1 for the groundstate. This localization length is 
slightly larger than the value predicted from the asymp- 
totic solution (£q — 1.23) of the simplified set of BdG 
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FIG. 8: (color online) Wavefunctions u for the zero mode and 
a nearly degenerate pair of the n — 1 modes at m = 1/2. 
Upper panel: u in logarithmic scale, showing the exponential 
localization of these modes, which is found to be weaker for 
the excited states. Note the hybridization of the zero-modes, 
while each of the excited states is localized near a single vor- 
tex. Lower panel: showing the Bogoliubov functions u(r) and 
in dashed v(r) for the same modes. Note the phase shift be- 
tween u and v for the excited state. The inset, again using a 
logarithmic scale, shows that the phase-shift remains constant 
far from the vortex core. 



equations (34) that we discussed in section VA As the 



size of the vortex is reduced, and the chemical poten- 
tial is increased, the localization length converges to the 
asymptotic estimate. 

Secondly, another interesting feature of the n > 
modes is the occurrence of a phase shift between the u 
and v functions, which is displayed in the lower panel of 
Fig. [8j In the vortex-core, the two functions are oscillat- 
ing at slightly different wavelengths. Far from the vortex 
core, they have sinusoidal forms of the same decreasing 
amplitude, but with a constant phase shift. 



D. Experimental Consequences in Cold Atomic 

Gases 

To illustrate our findings in more familiar units, we 
devote this section to discussing the orders of magnitude 
of the different energy-scales involved in typical exper- 
iments for the known Feshbach resonance in cold 40 K 
gases. As shown in Appendix [XJ the scattering parame- 
ters which ultimately determine the physics of the p-wave 
superfluid can be evaluated explicitly from experimental 
data. This analysis yields the magnitudes of the cou- 
pling constant g of the two-channel model (11), as well 
as the dimensionless constant ci which sets the maxi- 
mum density of bosons, and thus the coherence length 
£ = \J 02/(1 + C2) in our dimensionless units. The scal- 
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FIG. 9: (color online) Using the values thus obtained, of C2 ~ 
14.4 and g « 1.4 • 10~ 46 Jmi, the subgap spectrum can be 
given as a function of the detuning S, resulting in Fig. [9] 
Spectrum of the subgap states, in units of Fermi energy, as 
a function of the detuning 5 for the Feshbach resonance in 
a 40 K gas. The graph includes subgap states with angular 
momentum \m\ < 10; higher angular momentum states would 
lie between the displayed subgap states and the bulk gap. The 
largest mini-gap is obtained at a detuning S ~ 0.062, and 
amounts to roughly 1% of the Fermi energy, as highlighted in 
the inset. 



which amounts to 68.1 for the specific case of 40 K, further 
assuming a confinement length of £ — 500nm. In partic- 
ular, the overall scale of energies becomes E = S~ 1 EeF, 
in terms of our dimensionless energies E. As S is only 
moderately large for 40 K, the mini-gap can be as large 
as 1% of the Fermi-energy. The maximum of ei is rather 
shallow, so this magnitude of the mini-gap is realised over 
a significant interval of detunings 5 f=s 0.05 . . . Q.'Sep- At 
larger detuning, the energy scale of the subgap states de- 
creases linearly with 6 and becomes exponentially small 
as 8 — > 2ep. To conclude with an example, for a trapped 
potassium gas with ep — 10kHz we predict that the mini- 
gap can be of the order of 100Hz at a detuning frequency 
of 8 1500Hz. 

This energy scale for the mini-gap is still rather small, 
and we now consider mechanisms to further increase t\. 
As energies are scaled with 5~ 4 , we can read off from 
Eq. (30) how this can be achieved. In particular, de- 



ing of the axes relative to the dimensionless units used 



in the bulk of this paper is set by the factor S [see (30)], 



creasing the perpendicular confinement length of the gas 
£ represents a simple means to enhance ei slightly. Fur- 
ther improvements can only be achieved by using a dif- 
ferent Feshbach resonance with larger coupling constant 
g and smaller C2. However, as very small C2 will result in 
large vortex cores, this parameter should not be smaller 
than about 1. To visualise the scaling of energies as the 
vortex size £ is varied, Fig. 10 indicates the dependency 
the number of bound states below the bulk gap depends 
on the dimensionless parameters of the problem. Besides 
increasing the overall energy scale, reducing S also shifts 
the maximum mini-gap towards larger values of the de- 
tuning, which would be easier to stabilise experimentally. 

In cold atomic gases, the presence of the subgap states 
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FIG. 10: (color online) Summary of how the number of 
Caroli-de Gennes-Matricon subgap states varies with the size 
of the vortex cores £ and the chemical potential p, [in dimen- 
sionless units according to Eq. . Typical type II supercon- 
ductors yield £ ~ 1, while f behaves according to Eq. (29 1 in 
cold atomic gases. 



can be probed experimentally by using RF spectroscopy, 
in the same way that had been proposed as a probe for 
the zero modesP^ The idea is to probe the amplitude for 
resonant absorption of an RF photon leading to a tran- 
sition to a different hyperfine state, and therefore pro- 
jecting the affected atom out of the spin-polarized Fermi 
gas. The bulk signal amounts to an absorption edge at 
hu> = -Ehyporfinc — A* + ^-B , and the presence of the zero 
modes in vortex cores results amounts to the addition 
of a series of linearly spaced peaks starting at the lower 
energy hu> = £%perfine ~ A< whose intensity decays on 
approaching the absorption edgeP^ The Majorana mode 
yields a sequence of peaks by coupling to different states 
of the continuum. For trapped gases, these are spaced 
by the trap-frequency ujj_. Finite energy subgap states 
result in additional series of absorption peaks occurring. 
In order for these different signals to be well separated, 
it is required that the spacing of peaks in each series be 
small compared to the typical energy gap between sub- 
gap states, i.e., that the trap frequency be smaller than 
the typical spacing of the subgap states. 



VI. CONCLUSIONS 



As the number of subgap states scales roughly as the 
ratio of the bulk gap to the chemical potential, we argue 
that the regime of small chemical potential, or strong 
coupling should be most suitable to overcome these issues 
and maximize the mini-gap to the first subgap state. 

We note in passing that unconventional realizations 
of p-wave superconducting order where the magnitude 
of the order parameter is set externalljI^Hlo] ma y re q U j re 
different strategies to optimize the topological protection 
with regard to the subgap states. 

In particular, we study the case of atomic Fermi gases, 
where this regime can be reached when tuning the system 
close to a Feshbach resonance. We solve the Bogoliubov 
de Gennes equations for strong interactions in narrow 
Feshbach resonances, and calculate how the relevant sys- 
tem parameters, namely the chemical potential and the 
coherence length, depend on the detuning. Unlike typical 
type II superconductors in the BCS regime, the size of 
vortex cores can be tuned in cold atomic gases, allowing 
large vortices. 

As Kopnin and Salomaa's perturbative solution of the 
vortex states is not valid in small chemical potential, we 
confirm numerically that the energy of the subgap states 
remains finite in small chemical potential; the subgap 
states are ultimately absorbed into the continuum above 
the bulk gap at chemical potentials of /i ~ 0.5 roA^. 
Among the most interesting aspects of the structure of 
the vortex states we considered, we note a small split- 
ting between the states associated with vortices and an- 
tivortices. In large vortex cores, we characterize the na- 
ture of additional subgap states with a non-zero radial 
quantum number. Unlike the case of s-wave supercon- 
ductors, these branches are not symmetric under rever- 
sal of angular momentum. In the particular case of 40 K 
gases, we find that the optimal value of the detuning 
5 w 0.124e^ yields a maximal subgap energy of the order 
of O.OIcf- Mechanisms for increasing this energy scale 
were discussed. 
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0449521 and PHY-0904017 (VG). 



To summarize the main results of this paper, we have 
shown that the presence of subgap states complicates 
the use of p-wave superfluids as the medium underlying 
a topological quantum computer. While topologically 
protected operations are still possible, as information re- 
mains localized at a single vortex core even as transitions 
involving the Majorana zero mode and finite energy sub- 
gap states occur, the state of a qubit depends additionally 
on the parity of the number of excitations of the subgap 
states. This disqualifies some known strategies for the 
read out of such quantum bits, especially for neutral su- 
perfluids which are the main focus of this paper. 



Appendix A: p-wave Feshbach resonance parameters 



In this appendix, we derive the parameters of the p- 
wave Feshbach resonance in 40 K. To do this, we need to 
restore all physical units, thus unlike in the rest of this 
paper where we set K = 1, here we restore H in every 
formula. 

P-wave superfluids with a Feshbach resonance are 
characterized by two parameters,^ the coupling constant 
g and the ultraviolet cut-off momentum A. Instead of A, 
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it is more convenient to talk in terms of the dimensionless 
combination 



(-2 



% 2 _ 

3tt 2 



mr g 2 A 



/J 4 



(Al) 



In the context of the current work, c 2 sets the coherence 
length of the superfluid, while g sets the scale of the gap 
function Ao. Here, we briefly review an argument from 
Ref. [43] describing how c 2 and g can be extracted from 
the detuning between closed channel bosons and open 
channel fermions in a two-channel model. We then apply 
this method to the experimental data from Ref. [57] to 
extract c 2 and g for 40 K. 

The relevant physics is that of elastic scattering of two 
atoms, as described by the p-wave scattering amplitude 
3/i(fc)Pi(cos0), with the wave vector dependency 



fi(k) 



— - + ck 2 — ik 3 



(A2) 



Crucially, the scattering volume v and the prefactor of 
the second order term c can be extracted from experi- 
ments and from precise numerical modeling of the Fesh- 
bach resonance. On the other hand, we can relate v and 
c to the parameters g and o%. Indeed, we know thal!^ 



6irh 2 , . &Tih A 
- (e - const) , c = — (1 + c 2 ) . (A3) 



Trig 



m*g* 



In turn, the energy eg, which physically represents the 
Zeeman energy splitting between the open and closed 
channels can be quite generally written as 



(A4) 



where [1b is the Bohr magneton and a is a dimensionless 
parameter controlling the Zeeman splitting and which we 
take to be close toa«2 (its precise value depends on the 
physics of Fcshbach resonance) . This allows us to write 



dv- 1 



Qiraiigh 2 
mg 2 



6ttH , . ... 
■ -Wl + c 2 . A5 
m z g z 



Solving these for g and c 2 gives 



9 2 = 



67ra /i^/i 2 



f,2 dv- 1 
n dB 



(A6) 



In Ref. [57l c and are given as a function of B for 
the p-wave Feshbach resonance occurring at the magnetic 
field of 198.4 Gauss for the hyperfine state |/, m/) = 
1 9/2, —7/2) of potassium 40 K. Using their data, we can 
find the values of c and dv^ 1 /dB at the resonance and 
substitute them into ( A5 1. This gives 



c 2 « 14.4, g 1.40 • KT 46 Jm5 . 



(A7) 



We can use the values found here to see at what de- 
tuning n becomes of the order of mA§/?i 2 , the units of 



energy controlled by Ao- To do this, we rewrite (29 1 with 
h reintroduced [cf . Eq. ( 29 )] 



1< 



2tt(1 + c 2 ) 



m 2 g 2 1 



* s- x 



1 



(A8) 



where we again introduced the parameter S from ( 30 1 , 
with h inserted as needed, 



S = 27r(l + c 2 ) 



m 2 g 2 



(A9) 



and as before x = 5/(2ep). Comparison with (A3) gives 

ci 



S = 



(A10) 



Substituting the value for the relevant parameter from 
Ref. [37| (and using I ~ 500nm for the transverse width), 
we find that S « 68.1 and 



-^=68.1-^ 
mAjj 1 — x 



(All) 



We can estimate from here that in order to achieve the 
regime where \x is of the order of mA§//i 2 , i.e. roughly 
where the subgap states disappear, we need to keep 



x = 0.014. 

2e F 



(A12) 



A more accurate calculation of the proposed target value 
of the detuning is given in section |VD| above, based on 
the evaluation of the full subgap spectrum. 



Appendix B: The subgap states for fi ^> 1: Kopnin 
and Salomaa's approach 

For completeness, this appendix includes a pedagogi- 
cal summary of the solution to the BdG equations in the 
presence a single vortex |5]), as first proposed by Kop- 
nin and SalomaaP We specialize to the case of vortic- 
ity k — 1, and point out that the BdG equations ^ 
resulted from an Ansatz for a state with angular mo- 
mentum (L z ) = to around the axis of the vortex. In 
particular, we highlight the consequences of the initial 
assumptions as we proceed through the analytical solu- 
tion, with the assumptions being that 



(Bl) 



The solution proceeds in three steps. First, we note that 
for large p,, the coupling of the equations ^ for u and v 
is weak, and the solution is expected to be close to the 
decoupled equations given by the terms in square brack- 
ets only, which are solved by Bessel functions. Therefore, 
we proceed with the following Ansatz as the zeroth ap- 
proximation 

u(r) = H$(qr)f(r), v(r) = H^(qr)g(r), (B2) 
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where we introduce q = \/2ji ^> 1, and ffW are the Han- 
kel functions of type i € {1,2}. To simplify notations, 
we explicitly state the following calculation for the case 
i = l. For constant / and g we solve the parts of the 
equations in square brackets. More generally, we may al- 
low these functions to be slowly varying, on distances of 
the order of 1 (as opposed to the Hankel functions which 
vary over distances ~ l/g)- 



Substituting the Ansatz (B2 
and assuming that r 3> l/q, 



into the equations §5§ 
re Hankel functions can 



be replaced by their asymptotic expansions H(r) ~ 
e ±lr I \fr. In particular, their derivative is essentially pro- 
portional to the function itself times i, up to a small cor- 
rection, namely ^H(qr) = iqH — ^H. Wc find the fol- 
lowing equations (no approximations, except the asymp- 
totic expansion for the Hankel functions) 



-h 2 {iqf + f 



m f 
r 



hh'f 



~Q H q 

2 J 2r 



hh'g = Ef, 
1 

Eg. 



(B3) 



This new set of equations ( B3 ) can once more be an- 



alyzed in a perturbative approach. At the first level of 
approximation, we consider only the terms of order q (in 
particular, dropping E which is of the order of 1) 



iqf' + h 2 iqg = 0, 
-h 2 iqf + iqg' = 0. 



(B4) 



The solution to these equations yields a decaying be- 
haviour at r — > oo, which reads 



■ J-J-dr'foV) 



9 



-/J-rfr'ftV) 



(B5) 



Note that this approach does not capture the dependency 
of the localization length on the energy, as discussed 
above in section V A| To proceed with the perturbative 
solution of Eq. (B4), we study corrections at the next 
order in q, taking the following AnsatJ^"^ 



/ «L(r) + *^,,«-L(r) + *^ ! 



9 



(B6) 



where we introduced a shorthand notation for the local- 
ization factor 



L(r) = e -K dr ' h2{r,) . 



(B7) 



Substituting this last Ansatz into the equations ( B3 ) and 



equating terms which are (/-independent [the terms pro- 



portional to q cancel due to (B4|], one obtains 

E 



V4 - h ip2 



2 



+ = [-E- 



-h 2 



L{r) 
L{r). (B8) 



We look for a solution in the form 



ipx = a(r)L(r) + P{r)L~ l {r) 

4'2 = -a{r)L(r) + p{r)L-\r) (B9) 



The resulting equations for a and /? are a coupled system 
of differential equations 

a'L(r)+ p'L-^r) = (e - y + h 2 ™^j L{r) (BIO) 

a'L{r)-P'L-\r) = (-E - ^ - h 2 ™\ L(r), 

which can be decoupled by taking the sum and difference 
to yield 



' = -h 4 



P' - (E + h 2 (r)™)L 2 (r). 



(Bll) 



Integrating these equations, it follows that 

a = - [ dr'h 4 (r') (B12) 
Jo 

= - J™ dr' (E + h 2 (r') ^) e- 2 fo' dr"h 2 (r")^ 

As the radial profile of the gap function h(r) goes to 1 at 
large distances this solution implies that ip\ and ip2 are 
well behaved at large distances (i.e., they vanish). 

Finally, we need to construct a solution to the full 
equations, which is regular at the origin. This is achieved 
by matching boundary conditions on the set of solutions 
that we have generated. The imaginary part of the Han- 
kel functions which we used in the Ansatz ( B2 ) is diver- 
gent at the origin. A regular solution can be obtained 
by a superposition of the solutions expanded in 
and H^ 2 \ The solutions using are identical to the 
present discussion up to flipping the sign of i in ( |B3[ ). 
For the superposition of the corresponding solutions to 
be regular at r — > 0, the ensuing criterion is that /(0) 
and g(0) are the same for either of the two Hankel func- 
tions, such that the sum yields a Bessel function which 
is regular as r — > 0. 



Checking the Ansatz ( B6 1 , the condition is equivalent 
to ■0i(O) = "02(0) = 0, as both functions enter with a 
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factor of i in (B6|, so they differ by a sign between the 
two Hankel functions. We always have a(0) = 0, but in 
addition it is required also that j3(0) — 0, which yields 



dr \E + h 2 (r) 



0. 



(B13) 



This yields the harmonic spectrum for the subgap states, 
E m = — muJo, with the spacing spelled out in Eq. 
One characteristic of this solution is that the energy of 
the subgap states depends crucially on the behaviour of 
the gap function near the vortex core. Indeed, if we had 
h =const, then ujq goes to infinity, as the integral over 
r' in the numerator of (|9| diverges in the origin. The 
only exception to this rule is the zero mode at angular 
momentum m=0. For it, as can be seen directly from 



(B13), this divergence is not relevant, indicating the ex- 



ceptional topological protection of the zero-mode. 



Appendix C: Matrix elements between subgap states 

To make the nature of the matrix elements coupling 
the Majorana zero-modes and other subgap states more 
explicit, let us rewrite the single electron creation oper- 
ators <v in terms of the Bogoliubov transformed basis 
using the inverse of the unitary Bogoliubov transforma- 
tion 



o f (r) [«n(r)S„ + <(r)St 



N v /2 



J2 [^(r)f j+U f(r)f] 



(CI) 



for an even number of vertices Ny. (In the case of odd 
Ny an additional Majorana mode is found at the edge.) 
Here, we distinguish the contributions of the propagating 
Bogoliubov quasiparticles b n , the (pairs of) zero modes 
Tj and further subgap states c v<m localized around the 
vortex v and numbered also by their angular momentum 
to. To show that both the left-hand and the right-hand 
side of Eq. (CI) are fermionic, we have expressed the 



zero-modes in terms of fermionic operators obtained by 
pairing up individual localized Majorana zero-modes as 



r, 



2j -*72j-iJ, 



(C2) 



with the associated Bogoliubov functions u9(r) = 
u 2j,o( r ) + * u 2j-i,o( r ) an d likewise for v®. However, given 



the property of the zero modes, 
bution can be rewritten as 



N v /2 



V J,0 



their contri- 



[^(r)f,+ u °*(r)fj] =£^ (r) 7 . 



(C3) 



which does not explicitly show the fermionic nature of 
the sum. We now apply the more compact notation of 
Eq. |C3|. Using Eq. (ICll) to express a scalar potential 



V = J d 2 rV(r)a^ (r)a(r) yields a bilinear form involving 
terms mixing 6's, c's and 7's. The resulting expression 
includes, most prominently, the matrix elements between 
the states localized around the same vortices 



d2rV ( r ) 



(C4) 



t ( r )cL 



(r)> + 



t ( r )4,m) 



As the zero-mode and subgap states all have support in 
a limited region around the vortex core, any disorder po- 
tential that varies on the scale of the vortex is bound to 
induce sizeable matrix elements for processes involving a 
single Majorana operator such as & v % or %c v . We have 
included explicitly only those processes which are cou- 
pling states at the same vortex; other matrix elements 
are small assuming well separated vortices due to the ex- 
ponentially decaying tails of the vortex states (c.f. Eq.[7|, 
as we discuss in more detail in the numerical part of this 
paper. 



Appendix D: Form of the BdG equations on the 
sphere 

We brief ly summarize the equations derived by Kraus 
et a/! 26 * 27 ! and transcribe them into the dimension- 
less units used in this paper. When expressing the 
Bogoliubov-de Gennes equations for p-wave superfluids, 
we choose a phase winding of — 2it for order parameter^ 
requires q = — |, corresponding to the basis functions 
in the presence of a single monopole fluxj 26 * 27 -^ the re- 
maining quantum numbers I, m span a complete basis of 
eigenfunctions of this problem with I = \q\ + s, \m\ < I 
(s — 0, 1, . . .). The BdG equations thus have the matrix 
form 



H(lm)(l'm<) 



E n 



( D1 ) 

where the bar distinguishes indices that relate to the 
function v and primed indices are summed over. We 
study the vortex pair consisting of a vortex at the north 
pole and its antivortex at the antipodep^with the vortex 
field 



R/^sin0 



[1 + (£/<<: sin <9) 2 ] 2 



o ir t> 



(D2) 



which conserves the angular momentum to = (L z ) as a 
good quantum number. The formulation in Ref. [57] also 
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introduces a pairing range £ p which we take to be much 
smaller than all other length scales. The off-diagonal 
(pairing) matrix elements of the resulting only couple the 
angular momenta m = M and fh = 1 — M, i.e., Eq. (Dl ) 



separates into a block-diagonal form indexed by M with 
remaining matrix indices I, I'. The entries arej^with Di 
as in Ref . |27| 



H, 



(lm) (I'm' 



n 



(lm) (I'm' 




(D3) 



(2J + 1)(2J + 1) 



[A + (-i)' +r A-']E^ 2L + 1 /i (° 4 ) 



I I' L 
1/2 -1/2 



/ V L 
-m m — 1 1 



where the last line denotes the 3-J symbols for the cou- 
pling of angular momenta, and fY are the expansion co- 
efficients of the vortex field Fy in regular spherical har- 
monics 



Generally, the radius of the sphere needs to be cho- 
sen much larger than the size of the vortex core. For 
our numerical calculations, we used a maximum cut-off 
( max = 400. To minimize finite size effects, it is desirable 
to choose the Fermi angular momentum lp as large as 
possible. Convergence of the numerical scheme depends 
crucially on the ratio of R/£, and the maximal radius and 
Fermi angular momentum were determined in each case 
by the requirement of convergence of the eigenvalues. 

For large L, it becomes difficult to calculate the expan- 



sion parameters ( D5 ) explicitly by numerical integration. 
Instead, we fit t 



le values of fi,[R, £] for L > 100 by a 



function of the form 



f^[R,Z]=a[R,Z] 



D -b[R,t]L l L c[R,i] 



(D6) 



i(Sl)*F v {0), 



(D5) which yields an excellent approximation. 
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